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Some  properties  of  matrix  sign  functions  derived 

from  continued  fractions  ,,  »-M 
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!  V 


Abstract:  This  paper  proposes  an  alternate  representation  of  a  matrix  sign  function  based  on  an  irrational 
function  described  by  a  continued  fraction.  The  properties  of  the  continued  fraction  and  the  truncated 
continued  fractions  are  investigated.  Also,  new  algorithms  for  computing  the  matrix  sign  function  are 
developed.  The  matrix  sign  function  is  then  extended  to  a  generalised  matrix  sign  function  for  directly 
solving  discrete-time  system  problems. 


1  Introduction 

Since  Roberts  [i  |  initially  introduced  the  matrix  sign  function 
and  its  applications  to  linear  systems,  many  applications  for 
solving  system  problems  have  been  developed  [1.  2.  3|.  A 
Newton- Raphson  type  algorithm  proposed  by  Roberts  ( 1  ]  and 
an  improved  algorithm  by  Balzer  [4]  have  been  used  as 
standard  algorithms  for  computing  the  matrix  sign  function. 

One  main  feature  of  the  matrix  sign  function  is  that  it 
preserves  the  eigenvectors  of  the  original  matrix.  This  property 
is  useful  both  for  studying  the  eigenstructures  of  matrices 
and  for  developing  applications  to  engineering  problems. 

The  matrix  sign  function  S  of  a  square  matrix  Ae C"xn 
with  Re  (o(A ))  =£  0  is  defined  by  [1 1 

.S'  A  sign  (4)  =  2  sign* (A  1  -  /„  (Id) 

where  /„  is  an  n  x  n  identity  matrix  and 


A  recursive  scheme  for  computing  the  matrix  sign  function, 
given  by  Roberts  1 1  ]  and  improved  by  Balzer  [4] ,  isasfollows: 


siin'(A)  =  ■-*—  <£  (A/„ -A)~‘dX 


Sg  +  i  =  akSk  +0  kSk' 


S0  —  A 


<;  is  a  simple  closed  contour  in  right  halfplane  of  X  and  encloses 
all  the  right-halfplane  eigenvalues  of  A . 

following  the  definition  of  eqn.  1 .  if  A  has  a  Jordan  form: 

7  -  block  diag  [7,.7.)  i  7,  *  7  (2a) 


A  =  MJM~  (2b) 

and 

.S'  =  sign  (A)  =  M  (sign  (7*)  *  sign  (7.)|  AT' 

=  M\ /„,  *  (  lni)\M-'  (2c) 

where  7,eCnU"'.  7.6Cn!l'"!.  and  n  =  ril  +  rt2  arc  the 
collection  of  Jordan  blocks  with  Re(a(A))>0  and  Rc(a(A)) 
<  0.  respectively.  A/eCnxn  is  a  modal  matrix  of  A . 

The  extended  matrix  sign  function  S  of  A  including  Re(o 
(A ))  =  0  is  defined  by  13] 

S  *  sign  (A)  =  Af  (sign  (7t)  ®  sign  (7.)  ©  sign  (70))  AT1 

=  A/(/„i  ©  (-/n2)  ®  0„j |  A/'1  (3) 

where  7neC"3x’,'1  is  the  collection  of  Jordan  blocks  with  Re 
(o(A))  =  0,  0„3  is  an  n3  x  n3  null  matrix,  and  n)  +  n 2  + 
«3  =  rt. 
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0th  +  0fe  =  1  and  Urn  ak  =  lint  0h  =  J . 

ft  -*  00  k~*  « 

if  Re  (o(A ))  )0  (4) 

The  algorithm  for  the  extended  matrix  sign  function  is  given 
by  |3) 

4o  =  Sh*  i  +  •Sfe  *  i  (5a) 

where  and  5^*,  denote  the  (k  +  l)th  iteration  of  the 

matrix  sign  algorithm  in  eqn.  4  using  So  =  A  +  el„  and  S„  = 
A  —  e/„ .  respectively,  e  is  given  by 


where  0  <p  <  1  and  A°  is  the  Drazin  inverse  of  A  [3] . 

The  algorithm  in  eqn.  4  is  known  to  be  a  Newton-Raphson 
type  and  is  often  used  as  a  standard  method  for  computing 
the  matrix  sign  function. 

In  this  paper,  we  define  an  alternate  representation  of 
matrix  sign  functions  based  on  an  irrational  function  described 
by  a  continued  fraction.  The  properties  of  the  irrational 
function  and  the  convergence  of  the  truncated  continued 
fractions  arc  investigated.  This  leads  to  new  algorithms  for 
computing  the  matrix  sign  function.  It  is  shown  that  the 
standard  algorithm  in  Eqn.  4  is  a  special  case  of  the  new 
algorithms  derived.  The  matrix  sign  function  is  then  extended 
to  a  generalised  matrix  sign  function  for  directly  solving 
discrete-time  system  problems. 

2  Scalar  sign  function  > 

In  order  to  develop  a  new  algorithm  for  computing  the  matrix 
sign  function,  we  define  a  new  representation  of  a  (scalar) 
sign  function  as  follows. 

Definition 

The  (scalar)  sign  function  of  a  complex  value  XSo(A)  is 
defined  by 

|+I  if  Re(X)  >  0  (ha) 

sign  ( X)  = 

l-l  if  Re(X)  <  0  (6 b) 

f^j  if  Rc(X)>0  (6 ic) 


if  Re(X)  <  0 


83  07  27  050 


g(\) 

X 


if  Re(X)  >  0 


g(X) 

X 


if  Re(X)  <  0 
where  g (X)  is  defined  by 

X  if  Re(X)  >  0 


(6e) 


(6/) 


gW  = 


l-Xif  Re(X)  <  0 


Note  that  Re(X)  =  0  is  not  included  in  the  definition. 

When  X  is  a  real  value,  a,  then  g(X)  in  eqns.  6 g  and  bh 
is  obviously  an  absolute  value  function  and 


g(o)  =  VP 


(7) 


In  a  similar  fashion,  when  X  is  a  complex  function,  g(\)  may 
be  defined  by 


£(X)  =  VP 


(8) 


with  proper  selection  of  branch  cut  to  match  the  definitions 
of  eqn.  7  and  eqn.  8.  To  derive  the  functiong(X)  in  eqn.  6  or 
eqn.  8,  we  consider  a  square-root  function  f:z  -*  VP that  has 
the  branch  cut  on  the  negative  real  axis  and  the  first  Riemann 
sheet  with  |arg(z)|  <  n  as  the  domain  of  z.  The  continued 
fraction  expansion  of  VP s  given  by  (5) 


/(z)  =  VP  =  1  +  — ' 
2  + 


z-  1 


(9) 


2  4 


z-  1 


In  order  to  study  the  domain  of  z  so  that  the  irrational 
function  /(z)  can  fully  be  described  by  the  continued  fraction 
expansion  in  Eqn.  9,  we  investigate  the  properties  of  the 
continued  fraction. 

Define  the  Xth  truncation  of /(z)  as /*  (z),  or 


/*t(z)  =  1  +2/;.,  (z) 

'where  /5  (z)  =  0  and 

*(*-») 


fk  (?)  ± 


1  + 


H^-i) 


i  + 


Hz-D 


Hz-D 

The  recursive  forms  for  a *  and  bk  are 
0*  =  «*-i  +  Hz- l)«fc-2  «-i  =  1. 


(10/7) 


00  =  o 
(110) 


bk=  b„.t  +Hz-i)*». 


b-i  =  0,  b0  =  1 


(11*) 


The  difference  equation  for  both  ah  and  bk  satisfies 
1  ~q~x  -  Hz-  1 ) <7 _ 2  =  0 


(He) 


where  q  ' 1  is  a  backward  shift  operator,  or  q  " 1  ak  =ak  .  , . 
The  zeros  of  the  equation  in  eqn.  1 1  c  become 


<7t  = 


1  +  \fz  ,  1  —  VP 

— - - and  q2  =  - - - 


(11  d) 


|X|e'*  if  -y  <  <t>  <  j 

(6 g) 

/k(z) 

JX if  J<<f><Y 

(6/t) 

From  eqn.  1 1  cl  we  observe  that 

|q,|>|q:|if-jr<arg(z)<7t  (He) 

The  general  form  for  eqn.  10a  (see  Reference  S)  is  as  follows: 

vH(l  +  VPf  4(1  -  n/z~)»| 

(i  +  VP )fc  -(i  -  VP )k 


kC2jz' 


v.r,  y 

. w  k  ^  i]  + 1  * 

J  =  0 


(12a) 


where  p  =  (A/2 1  and  r  =  ((X  -  I  )/2|  are  integers  for  lc  >  1 
and  j ,C2j  are  the  coetficicnts  of  a  binomial  expansion. 
Substituting  eqn.  1  !</  into  eqn.  1 2a  yields 


/*(*> 


Vz  Ui'i  +  q2  ) 

Hz)' 

‘A  -  </; 

1 

_(«£ \fc 

Wi  / 

=  \/z 


1  4  I 


k  >  1  and  </,  4  q2 

fk(:)  tor  k  =  1 . 4  arc  as  follows: 

/.(--)  =  I 

Z  4  1 

hi:)  =  ----- 


(10a) 


k>\ 


hi:)  = 

/*4  (-  )  = 


dr  -f 

z  4  3 
-  2 


4  6z  4  ) 


4z  4  4 


(1 2b) 

(13a) 

(13/r) 

(1 3c) 

( 1 3c/) 


Some  properties  of  the  continued  fractions  in  eqns.  9  and  10 
which  will  be  used  to  derive  the  g(X)  in  eqns.  t>  and  8  are  as 
follows: 


Property  I 

If  l<7 1 1  >  \q2 1.  then  we  have 
/(z)  =  lim  fk(z )  =  V? 


(14a) 


The  important  result  in  eqn.  14a  can  be  verified  by  a  ratio  test 
of  the  series,  which  can  be  obtained  by  expanding  [I  — 
(<7:/7t)*l  1  in  eqn.  12 b.  Also,  the  convergent  condition 
implies  that  the  domain  of  z  is  in  C*  where 
C  =  C  -  R  *  and  R"  is  the  negative  real  axis  R"  =  (— *>,  0] , 
or  that  the  complex  variable  z  must  be  -tr<arg(z)<». 
Thus,  the  function  f(z)  uniformly  converges  to  the  desired 


112 


IFhPROC.  lot  130.,  Pt  D.  No.  3  MAY  1983 


function  \Jz ,  and  \fz  can  fuUy  be  represented  by  the  con¬ 
tinued  fraction  if  z  €  C  *  . 

Property  2 

If  \Q\ I  =  k/j  I  and  z  =  0,  then 

f(z)  =  lim/fc(z)  =  0  (14A) 

k  — «■ 

Property  3 

If  l</i  I  =  I<7j  I,  |z|  #  0,  and  |arg  (z)l  =  jr,  then 

lim/fc(z)-°°  (14c) 

k  —  “ 

The  results  in  eqn.  14*  and  c  can  be  verified  as  follows: 

When  li/il  =  k?2l.  from  eqn.  1  \d  we  have  |l+\/z|  = 
II  —  y/z\.  If  z  =  0,  then  we  have  ah  =  — (X/2X1/2)*  and  A*  = 
(X  +  1X1/2)*  in  eqn.  11.  Thus,  /*(z)  in  eqn.  10a  becomes 
1  /(X  +  1 )  and 

/(- )  =  lim  /*(z)  =  lim  — 7—  -*-0  (14d) 


where  C  ^  C  -  land  I  is  the  entire  imaginary  axis.  From  the 
properties  shown  in  eqn.  14,  the  domaip  of  z  is  in  C "I there¬ 
fore  the  domain  of  X  must  be  in  C .  In  other  words,  the 
function  g(X)  converges  if 

X  =  |X|</*  and  |0|  *  ^  (ISA) 

Based  on  the  convergent  condition  derived  in  eqn.  14a,  wc 
have  the  desired  function  g(\)  in  eqn.  6  as  follows: 

*(X)  =  Jimg,,(X)  =  n/X1 

s/IXIV21*  =  =  X  if 

VixFe77®^  =  |X|e>(4’-ff)=  -X  (16a) 

,  rr  3rr 
if  -<</><  —  (16*) 

where 


The  function /(z)  converges  to  zero  if  2  =  0. 

On  the  other  hand,  when  z  is  a  nonzero  negative  real,  then 
\/z=jcj.  Thus  1  +>/z  =Vl  +  w2  e**  and  1  —  \/z  =Vl  +  w3 
where  0  =  tan' 1  to.  Therefore,  we  have 


I  +  *'J2**  w 
)h(2)  ~  ,U  1  -e-nk*  ~  tan(X0) 

The  function  fk(z)  diverges  as  z  €  (—  <*>,  0). 


(14c) 


Property  4 

The  poles  and  zeros  of  /„(z)  and  k  >  2  alternate  on  the 
negative  real  axis. 

From  eqn.  1 2  we  have  the  poles  of  /*(z) 


m 


1.2 


(X  —  1 )/ 2  for  odd  k 
(X  —  2)j2  for  even  k 


(14/) 


8k(k)  =  ft >(*)  with  z  =  X2  (16c) 

Thus,  the  scalar  sign  functions  defined  in  eqn.  6  can  be 
expressed  by 


s,gn(x)  s  i5> 


Slgn&>  M 


or 

sign(X)  ^  ^  =  lim  sign&^X) 
A  A?-* 00 


(17a) 


(17*) 


where 


signal*)  = 


X 

ft(X) 


(1  +X)*  -(1  -X)* 
(1  +X)*  +(1  -X)* 


signal  (X)  k 


ft(X)  _(1  +  X)*  +  (1  —  X)* 
X  (1  +  X)* -(1 -X)* 


(17c) 


(17t/) 


and 


and  the  zeros 


z 


m 


tan 


(m  -  $  )rr  2 
k 


*fc(X)  =  X 


(1  +  X)*  +  (1  -X)* 
(1  +  X)*  —(1  —  X)* 


k  =  1,2,..  . 


is  the  Xth  truncation  of  the  continued  fraction 


(17c) 


m  =  1,2,..., 


(X  —  1  )/2  for  oddX 
X/2  for  even  X 


(14g) 


Since  the  tangent  function  tan*  is  monotonically  increasing 
for  0  <  *  <  ir/2,  the  poles  and  zeros  of  fh(z)  for  X  >  2  alter¬ 
nate  on  the  negative  real  axis. 

Using  the  function  f{z)  and  the  properties  obtained  in 
eqn.  14  we  are  now  able  to  derive  the  desired  function g(X)  in 
eqn.  8  using  the  principal  square-root  function  f{z)  of  eqn.  9. 
Let  z  =  X2 ,  we  have 

fit)  ±  s(X)  =  Vx*.  xeC' 

X2  —  1 _ 

_.x2-l  (I5a) 

2  +  — rr — ; - 


g(X)  =  lim  gk(X) 

k~*°° 


=  1  + 


1*  - 


2  + 


2  + 


x2-i 


xeC'  (17/) 


3  Matrix  sign  function 

The  scalar  sign  function  derived  in  Section  2  can  be  extended 
to  a  matrix  sign  function.  For  this  extension  we  need  to 
investigate  a  matrix  function  generated  by  a  scalar  analytic 
function. 

Consider  a  matrix  A  €  C  "  x  B  with  spectra  a(A )  =  {X ( . 

X,}  where  l<n.  If  a  scalar  function  p(X)  is  analytic  at  X/. 
/  =  1 . I,  then  the  matrix  function  P{A )  generated  by 
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p(X)  can  be  defined  as  (6) 

P(A)  =  ~ 7  6  P(X)(X/„  -^)'1  dX  (18) 

2 m  Je 

where  c  is  a  simple  closed  contour  which  encloses  Xy,/  =  1 
The  matrix  function  described  in  eqn.  18  has  the  following 
properties  [6] : 

Lemma  1 

Let  A  be  defined  as  above  and  p(X),  <7(X)  and  r(\)  are  analytic 
at  X;-  6  o{A),  \j  =  1 . /<n,then 


Since  the  poles  and  zeros  of  /(z)  in  eqn.  1 4  alternate  on  the 
negative  real  axis,  the  poles  and  zeros  ofg(X)  in  eqn.  17 c  or 
in  eqn.  \ld  alternate  on  the  imaginary  axis.  As  a  result,  if 
o(A)€C\  or  no  Xy€o(/l)  exist  on  the  imaginary  axis,  then 
^(Xy),  /=  1,  .  .  .  ,  /  are  finite.  Thus,  from  lemma  2  and  the 
result  shown  in  eqn.  16,  we  have 


if 

if 


Re(Xh)>0 

Re(X*)<0 


(21a) 

(216) 


Since 


(i)  if  p(X)  =  k.  then  P{A)  =  kIn 

(ii)  if p(X)  =  X,  then P(A)  =  A 

(iii)  if  p(X)  =  q(X)  +  r(X),  then  P(A)  =  q(A)  +  r{A) 

(iv)  if  p(X)  =  q(X)  KX),  then  P(A )  =  q(A  )  r(A )  =  HA )  q(A ) 

(v)  if  p(X)  =  r(q(X))  and  r(\)  is  analytic  at  q(X,),  and  there 
exists  Xy  £  a(A ),  then  P(A )  =  r(q(A )). 


When  A  has  kj  Jordan  chains  of  length  corresponding  to 
Xy  e  a(A ),  /  =  1 . 1  and  X  tjh  =  m,  where  m,  is  the  multi¬ 

plicity  of  Xj,  then  A  can  be  represented  by 
A  =  MJM~l  (l9fl) 


where 

J  —  J  x  ®  Jj  ®  .  . .  ®  J\ 


(196) 


and 

J.  =  /y,  ®  Jjj  ®  . 


h*  = 


Xj  1 
0  Xy 
0  0 


0 

1 


[_0  0  0 


Jlhj 


0 

0 


ec  >jh x  ijh 


L 


X/ 


(19c) 


(19  d) 


i'i'Xk)  = 

we  have 

g(Jkh)  = 


f  1  if  Re(Xfe )  >  0 
l-l  if  R  e(Xk)<0 


(21c) 
(21  d) 


Jkh  if  Rc(X*)>0  (21e) 

\-Jkh  if  Re(X* )  <  0  (21  f) 

Using  the  result  of  eqns.  21c  and  21/,  eqn.  20a  becomes 

g(A)  =  M[Jt  ®/2  ®  ...  *(-Jntl) 

©...  ®  (-/,)!  AT1  (21g) 

where  Re(Xk)>0  for  \<k<m  and  Re(X*)<0  for  m<k 
<  l.  Finally,  the  desired  matrix  sign  functions  can  be  obtained 
from  eqns.  21#  and  17  as  follows: 


sign(,,(4)  £  ^7  Xg(X)  '  ^  ^ 

=  A\g(A)Vx  =  M [/m j  ®  /mj  e 


W'  (22a) 


or 


sign(2)(/t)  ^  — <j)  ^ 


h  =  1 . kj  and  /  =  1 ,  .  .  .  ,  / 

Using  lemma  1,  a  matrix  function  g(A )  generated  by  g(X) 
becomes 


g(A)  = 

M ~ ;  <)>  g(X)  (X/ „  — 

J)-'  dXM~'  = 

MgfJW1 

= 

M {gfJjl  )  ®  g(Jt 2  )  ® 

©  giJjkj))M 

/  =  1 . 1 

(20a) 

where 

t-0 

h  ~  1 . 

kj. 

/•  =  1 . 1 

(206) 

gw>(X)  is  the  rth  derivative  of  g(X)  for  t  >  1 ,  and  £<0>(X)  =  g(X). 
Ht  .  a  shift  operator,  is  the  r*  x  r*  matrix  with  null  diagonal 
entries  and  all  Is  on  the  super  diagonal  entries.  Since #(/)  is  an 
upper  triangular  matrix,  we  have  the  following  result  [6] . 

Lemma  2  ,  ,  „ .  . 

Given  A€Cn*n,  <j(A)={\jJ=  1 - ,/</»}  and  g(X)  is 

analytic  at  each  Xy  €  <r(A ),  and  g(Xy)  is  finite,  then  oig(A ))  = 
&Xy),/=  1 . /<«}. 


=  A''  |g(4)l  (226) 

and 

sign(A)  =  sign(,)(/l)  =  sign(:>(A)  (22c) 

Thus  using  lemma  1  and  the  results  in  eqn.  17.  we  have  the 
following  theorem. 

Theorem  I  _  „ 

Given  a  matrix  AeCnXn  and  o(/l)6C.  the  matrix  sign 
function  of  A  can  be  described  by  a  matrix  continued  fraction 
[7 J  as 

sign(/4)  =  A[ln+(A2  -/„)  [2 /„  +  {A1  -  ln)  x 

[2I„  +  (4 2  -/„)[...]■'  ] _I  ]"'  l'1  (23a) 

=  A-'[yn  +  (A2  -In)[21n+(A2  -/„)* 

[21n  +  (A>  -ln)\... (236) 

Corollary  1 

The  matrix  sign  function  of  A  as  defined  in  Theorem  1  can  be 
approximated  by 

sign^ (/4 )  =  [(/„  +  A  )"-(/„-  A  )"  ]  (( /„  +  A  )k 
+  {ln-A)k}-' 


114 


/WWW..  I  ol  HO.  ft  n.  So.  3.  MAY  1983 


=  LIJ*Cjy+1/»^«  X  kC2iA» 


k  =  1,2,... 


signg>>(4)  =  [(/„  +  A?+(ln-A)k\  [(/„  +4)* 

-(In- A?]'' 

[W  i[M  r 

-  X  kCi}A2i]  [  £  fcCjy+  \A2>* 1 

J  -  A  J  1  i  _  A 


k  =  1,2,...  (246) 

For  convenience  of  application  we  list  some  matrix  sign 
functions,  sigr$J  (A)  for  k  =  1 . 7,  as  follows: 

sign  $(4)  =  A~'  (25  a) 

sign}2/,  (4)  =  (A2  +  /„)(24)''  (25  6) 

stgn}§}(4)  =  (34 2  +  /„)(43  +  34)-'  (25c) 

sign^{ (4)  =  (/l4  +  64 2  +  /„ )  (4A 3  +4/1)-'  (25c/) 

sign{/,(4)  =  (54  4  +  \0A 2  +  /„)x 

(4  s  +  104  3  +  SAY'  (25e) 

sign}/,  (4)  =  (-4 6  +  15.4 4  +  15/1 5  +/n)x 

(bA 5  +  204 3  +  6/1)-'  (25 f) 

sign}  V)  (A )  =  (7A6  +  3544  +  21A2  +/n)x 

(A2  +  21.4  s  +  35/4 3  +  74)“‘  (25*) 


Observe  that  sign}2}  (/l)  in  eqn.  256  is  the  standard  matrix 
sign  algorithm  in  eqn.  4  for  k  =  0  and  a*  =  0*  =  $ . 

For  a  system  which  contains  Re(o(4 ))  =  0.  the  matrix  A 
is  modified  as  shown  in  eqn.  5.  Then  the  proposed  matrix 
sign  algorithms  in  eqn.  24  can  be  applied  to  determine  the 
extended  matrix  sign  function  [3  J . 


4  Computational  considerations 

For  online  computation,  the  matrix  sign  function  in  eqn  24  or 
eqn.  25  can  be  used  to  approximate  the  matrix  sign  function 
in  eqn.  23.  However,  if  the  matrix  A  contains  both  large  and 
small  eigenvalues  in  modules,  for  a  large  k  of  sign}},)  (A)  or 
sign}*}  {A )  in  eqn.  24.  numerical  difficulty  might  occur.  To 
resolve  this  difficulty,  the  following  recursive  algorithms  are 
derived. 

From  eqn.  I  la  we  have 


sign}^,  (X)  = - 7  — r* 

1+V+V 

Rearranging  eqn.  26a  gives 

( 1  ~  X\*  _  1  ~  sign}j,))(X) 

\1  +xj  l+sign}/,(X) 


Rearranging  eqn.  26a  gives 


Let  k  =  k,k3.  we  have 


signW,  *,>(*)  = 


W 
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=  sign}!?,)  l sign}**) (X)]  (27a) 

=  sign}*,)  [sign}1**,  (X)]  (276) 

Letting  k  =  ktk2  ■  ■  .  km  and  using  eqn.  27a  repeatedly  yields 

signal  (X)  =  sign}**)  (sign}'*3,  [. .  |sign}/m,  (X)]  \ )  (27c) 

Similar  recursive  algorithm  can  be  derived  for  sign}*}  (X)  in 
eqn.  176. 

Theorem  2 

Recursive  algorithms  for  computing  the  matrix  sign  functions 
of  A  for  a(4)G  C  D  are 

s‘gn('4+  i> (A )  =  sign}}*,  [sign}},*)  (A )] ; 

s*gn(i)  (A)  =  A  (28 a) 

or 

sign}/*  + ,)  (4 )  =  sign (sign}/*,  (A )) ; 

sign},2}  (4)  =  A' 1  (286) 

where  nk,  ,  =/*rt*  for  k  =  1,  2 . /*  >  1  and  n,  >  1 

Remark  1 

Using  the  following  property  of  a  matrix  sign  function  sign(4) 
=  sign(4  “ 1 ),  we  can  set  the  initial  condition  of  eqn.  286  to  be 

sign}2/,  (4)  =  4  (28c) 

Remark  2 

The  standard  matrix  sign  algorithm  in  eqn.  4  is  in  fact  a 
special  case  of  the  matrix  sign  algorithm  derived  herein  by 
choosing  /*  =  2  in  eqns.  28 a  or  6.  For  example,  from  eqn. 
256  we  have 

sign}!}  (4 )  =  (4 2  +  /„ )  (24  )’ 1  =i(4+4‘')  (28c/) 

and 

sign}2/, (4  )  =  sign}2/,  (sign}2/,  (4 )] 

=  (44  +  64 2  +  /„ )  (44  3  +  44  )~ 1  (28c) 

The  result  in  eqn.  28c  is  identical  to  that  of  the  recursive 
algorithm  in  eqn.  4  using  k  =  1  and  a*  =  0*  =  1/2.  Other 
new  recursive  algorithms,  sign}/,  (4)  for  prime  number  k>  3 
in  eqn.  25  can  be  considered  as  basic  recursive  algorithms  fot 
computing  the  matrix  sign  function. 

Remark  3 

From  eqn.  26a  and  lemma  2,  if  X,  =  1  +  pei't>*.  0  <  |p|  «  1 . 
is  an  eigenvalue  of  4.  then  X(*  is  the  corresponding  eigen¬ 
value  of  sign}/,  (4 ): 


/  p  \k 

-  1  -2  "  eik*i 


(29a) 

US 


Then  we  have 


Similarly ,  if  \  =  —  l  +  pe^>  0  <  |p|  «  l ,  then  we  have 


Therefore,  the  order  of  convergence  in  the  neighbourhood  of 
±1  is  k.  Furthermore,  if  X(  =  ±  1 ,  then  the  corresponding 
eigenvalue  of  sign^J^  (A )  stays  at  ±1.  The  same  remarks  are 
true  for  sign{j^  (A ). 

Remark  4 

If  sign^)  (A )  or  sign^)  (A )  is  a  satisfactory  approximation  of 
sign  (4),  and  if  the  recursive  algorithm  in  eqn.  28  with 
constant  order  /*  and  number  of  iterations  m  is  used,  then  the 
number  of  iterations  needed  for  large/small  eigenvalues  is 
log,fcyV.  This  result  can  be  verified  as  follows: 

Let  N  <(fk)m ,  then  m  -  log^.V  (29c) 

Thus,  the  result  obtained  by  using  the  approximate  model  in 
eqn.  24  with  k  =  (fk  )m  —  N  is  equivalent  to  the  result  using 
the  recursive  algorithm  in  eqn.  28  which  has  a  constant  order 
fk  and  the  number  of  iterations  (m)  shown  in  eqn.  29c. 


5  Generalisation  of  matrix  sign  functions 

The  matrix  sign  function  defined  in  Section  3  can  be  viewed  ,s 
a  nonlinear  mapping  which  maps  the  eigenvalues  at  the  right 
and  left-hand  side  of  the  imaginary  axis  to  +  1  and  -  1 , 
respectively.  Also,  the  matrix  sign  function  preserves  the 
eigenvectors  of  the  original  matrix.  The  above  properties  are 
useful  for  examining  the  eigenstructure  of  a  matrix  and  for 
solution  of  control  system  problems. 

The  matrix  sign  function  can  be  generalised  to  map  the 
eigenvalues  of  a  matrix  to  +  1  for  those  eigenvalues  located  at 
one  side  of  a  simple  closed  curve  in  €  and  to  —  1  for  those  at 
the  other  side  of  the  simple  closed  curve. 

Theorem  3 

Let  l.C  C  be  a  simple  closed  curve  which  can  be  mapped  onto 
the  imaginary  axis  by  a  conformal  mapping  h(\).  Assume  that 

a  matrix  36C"X"  with  o(A  )  =  {X,.  i  =  1 . /}  exists  such 

that  o(A )  o  L  -  <t>  and  h( X)  is  analytic  at  X(.  Define 


5(A)  =  sign{h(A)]  *  Af  |/m,  ®  . . .  ®  l„m  ® 

(-/«„.,)  ®  (30 d) 

if  the  Jordan  decompostion  of  A  is  given  by 

A  =  M[Ji  ®  ...  ®  Jm  ©  Jmtl  ©  ...  ©  7/| AT1 

(30e) 

where  J(  is  a  generalised  Jordan  block,  and 
a(Jt)  =  {X, :  Re(/i(X,))  >  0  for  1  <  /  <  m  and 

Re(/i(X,))<0  for  m<i<l\  (30/) 

Proof 

Since  L  is  a  simp!,  closed  curve  in  C  and  h(k)  is  conformal, 
which  maps  L  onto  /co-axis,  the  whole  complex  plane  is 
separated  into  two  regions  Ci  and  C2  by  L,  such  that 
Re(Jr(X)) >0  for  X6C,  and  Re()r(X)) <0  for  X 6 € 2 .  Since 
/i(X)  is  defined  to  be  analytic  at  X,  €  o(A ),  i  =  1 ,  . .  . ,  /,  from 
eqn.  18  we  have 

h(A)  =  '  /i(X)(X7n -A)"1  <A  (31a) 

2  m  -Tc 


Assuming  that  o(A  )  1  /.  -<P  yields.  /i(X,)£  o(/r(A  )).  /=  1 . /, 

which  are  not  on  the  imaginary  axis,  or  Ji(Xj)  #/co,  u£R. 
Thus,  we  conclude  that  /r ( X, )  is  in  the  domain  of  sign  (X). 

Define  5(X)  -  stgn(/r(X))  (316) 

where 


sign(X)  =  — 


X 

jrtX)- 


(31  c) 


Thus 


S(  1 1  -  —  0  5<X|(X/n  .4)-’  J\ 

2m  -c 


1  :  g(h(X)l 

© 

2m -c  6(A) 


<X/n  -  AY'  d\ 


i  j  ink) 

2 m®  ,g(6(A)) 


(kln 


-  A  f 1  JX 


(317) 


5(A)  i 


1  r  gjh(k)) 
2m  ?  6(A) 


(X/n -A)'1  Jk 


or 


(3Ca) 


If  A  can  be  decomposed  into  the  form  of  "qn.  30e.  from 
eqn.  20  we  have 

h(A)  -  M\h(Jx  1  ©  ...  ©  h(Jm)  ® 

h(Jm.,)  *  .  .  .  ©  h(J,)\M~'  (3 le) 


sw  -  M 


h(k) 

mx)) 


(kl„  -A)'1  7X 


where 


S(X)  =  I  + 


XJ  -  I 


2  + 


XJ  -1 


2  + 


X* -i 


A[(l  +A)»+(1  -X)*] 
(1  +  X)*-(1  -xf 


(30 b)  Furthermore,  since  Re(/i(X,))>  0.  I  and  Re(/i(X())< 

0.  m  <i  <  I,  from  eqn.  22.  we  have 

5(A)  =  sign(/i(A)>  =  M|/m,  ©  ...  ©  l„m  © 

e  (“W1W-1  (31/) 

In  a  manner  similar  to  that  of  theorem  2,  we  have  the 
following  computational  algorithm  for  the  generalised  matrix 
sign  function. 

Corollary  2 

Let  h(X)  and  A  be  defined  as  in  theorem  3.  Then  we  have 
(30c)  5(A )  =  sign(6(/ ))  =  lim 5„.  (A)  (32«) 

lr_«,  " 
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where 

S«k.M)  =  Sfk  [5„fc(/l)l;  Sni (A )  =  6(A)  (32*) 

=  fk  nk  for  A:  =  1,2 . fk>l,n,>l 

and 

$*(*)  =  [(/„  +  h(A))nk  -(/„  -6(A))"*-]  x 

[('„  +  h(A  )Tk  +  Un  ~  h(A  ))"*  ] ' 1  (32c) 

for  k  =  1,2,... 

With  theorem  3,  we  can  develop  appropriate  matrix  sign 
functions  for  several  applications.  For  example,  in  discrete- 
time  control  system  problems,  we  usually  need  to  separate 
the  eigenvalues  of  a  matrix  A  by  the  unit  circle,  |X|  =  1, 
which  can  be  mapped  onto  the  imaginary  axis  by  a  class 
of  conformal  mappings  [8],  6(X)  =  (a\  +  6)/(cX  +  d)  with 
a  =  1 .  6  =  —  1 ,  c  =  1  and  d  =  1 ,  or 

6(X)  =  (33c) 

Replacing  X  in  eqn.  1 7 a  by  6(X)  in  eqn.  33a  gives 

S$,(X)  =  sign&(6(X))  =^i1  (33*) 

Also,  similarly,  eqn.  176  yields 

=  j _  |  (33c) 

Thus,  the  corresponding  generalised  matrix  sign  functions  for 
the  mapping  shown  in  eqn.  33c  become 

S\i\(A)  =  (Ak -/n)(Ak  +inyl 

k  =  1,2, . . .  (34c) 

or 

■saU-o  =  o*  +in)(Ak-inrx 

k  =  1,2....  (346) 

and 

S(A )  ~  lint  .^(A)  =  lim  S$(A)  (34c) 

k  —  00 

The  algorithms  derived  in  eqn.  34  for  computing  generalised 
matrix  sign  functions  can  be  directly  applied  to  solve  the 
algebriac  discrete  Riccati  equation  (see  Appendix),  discrete 
regulator  problem  (9) ,  and  discrete-time  stability  problem 
[10]. 

6  Conclusions 

This  paper  has  proposed  an  alternate  representation  of  the 
matrix  sign  function  based  on  an  irrational  function  described 
by  a  continued  fraction.  It  has  been  shown  that  an  irrational 
function  \/z  of  a  complex  variable  z  can  be  fully  represented 
by  a  continued  fraction  if  z  is  not  a  negative  real  value.  Also, 
the  poles  and  zeros  of  the  truncated  continued  fraction  alter¬ 
nate  on  the  negative  real  axis.  The  principal  square-root 
function,  f\z)  =  \fz,  is  then  extended  to  generate  a  matrix 
sign  function.  It  has  been  shown  that,  when  the  matrix  of 
interest  has  no  eigenvalues  on  the  imaginary  axis,  the  matrix 
sign  function  can  be  fully  described  by  a  matrix  continued 
fraction. 

Based  on  the  structure  properties  of  continued  fractions, 
new  recursive  algorithms  for  computing  the  matrix  sign 


function  have  been  derived.  It  has  been  shown  that  the  most 
commonly  used  matrix  sign  algorithm  is  a  special  case  of  the 
recursive  algorithms  derived  in  this  paper.  Finally,  the  matrix 
sign  function  is  extended  to  a  generalised  matrix  sign  function 
such  that  the  newly  developed  matrix  sign  algorithms  can  be 
directly  applied  for  solution  of  discrete-time  system  problems. 

For  continuous-time  systems,  the  proposed  matrix  sign 
functions  can  be  applied  to  solve  stability  problems  [3] , 
Riccati-type  and  spectral  factorisation  problems  [2] . 
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9  Appendix 

To  show  the  procedure  of  using  the  matrix  sign  functions  for 
solving  discrete-time  Riccati  equation  and  discrete  regulator 
problems,  we  consider  the  following  discrete-time  system: 

X(k  +  1 )  =  AX(k)  +  Bu(k)  (35a) 

y(k)  =  CX(k)  (356) 

where  A  eC"x  ",  SeC"*  m.  CeC**  "  and  A  is  nonsingu¬ 
lar.  The  infinite-time  performance  index  to  be  minimised  is 
given  by 

J  =  t  [XT(i)QX(i)  +  uT(i)Ru(i)]  (36) 

i"0 

where  Q  =  is  a  symetric  nonnegative  matrix  and  R  is  a 
symmetric  positive  definite  matrix.  Assume  that  the  pair 
\A ,  Q)  is  detectable  and  {A .  B)  is  stablisable.  then  the  steady- 
state  feedback  optimal  control  law  [8]  becomes 

u(k)  =-(/?  +  BTPByi  BTPAX(k)  (31a) 
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where  the  nonnegative  definite  matrix  P  is  the  solution  of  the 
following  algebraic  nonlinear  discrete-time  Riccati  equation: 


P  =  Q+AtPA-AtPB(R+  BtPB)'1BtPA 


(37  *) 


The  procedures  to  determine  the  optimal  control  law  using  the 
matrix  sign  functions  are  described  as  follows: 

Define  a  2n  x  2n  Hamiltonian  matrix  [8] : 


G  = 


A A'lBR-'BT 
QA~l  At  +  QA-lBR-'BT\ 

The  modal  matrix  of  G  and  its  inversion  are  defined  as 


(38a) 


M  ± 


Mxi 

AT1  £ 

Mn 

Ml2 

M2i 

M22 

Mj  i 

M» 

Mue  C"Xb;  /,/  =1,2 
MfjeCnXn-,  /,/  =1.2 


(38*) 


Thus, 


M~'GM  =  block  diag[A,  A'r]  A6C 


n  X  rt 


(38c) 


where  A(A~T)  is  the  Jordan  block  corresponding  to  the 
eigenvalues  of  G  outside  (inside)  the  unit  circle.  The  Riccati 
matrix  gain  P  in  eqn.  37*  can  be  determined  [8]  by 


P  =  MnMu' 


(39) 


P  in  eqn.  39  can  be  indirectly  computed  via  the  matrix  sign 
functions  as  follows: 


sign(G)  =  M  block  diag[sign  (A),  sign  (A  T)]  -M~ 1 

=  Mblockdiag[/n,-/B]Ar'  =  MTM~x  (40) 
where  sign(A)  =  /„ ;  sign(A~T)  =  —  /„;  /  =  block  diag(/„.~ 


/„]  and  /„  is  an  w  x  n  identity  matrix. 
Define  a  new  matrix  W ,  or 

W  |  sign(G)  +  /  =  M\TM~x  +3T1/] 


=  M- block  diag[2Afn,  —  2Mi2\ 


=  2 


MnMu 

—  MuMn 

A 

ft'u  Wn 

MnMn 

-MnMn 

IVji  W22 

=1,2 


Thus  P  in  eqn.  39  can  be  indirectly  determined  using  the 
partitioning  matrices  k^.and  Ifn  as 


P  =  =  (2M2lMn)(2MuMu)-' 


r2l  "'ll 

=  M21ATnl 


(42) 


As  a  result,  the  optimal  control  law  in  eqn.  37a  can  be 
obtained. 

In  practice,  an  approximate  sign(G)  is  often  used  to  deter¬ 
mine  the  P  in  eqn.  42.  The  matrix  sign  algorithms  for  com¬ 
puting  the  approximate  sign  (G)  (defined  as  Gy)  are 

sign(G)  =  _lim(G>-/Jn)(G^  +  /jBr' 


=  lim(G'  +  /2n)(G' 

j  —>oo 


f2„r‘ 


(43a) 

(43*) 

(43c) 


—  Gj  for  a  finite  j 
The  index  j  of  Gj  in  eqn.  43c  can  be  determined  when 

It  race  [(Gy)2 )  -  2m|/2m  <  e,  (43d) 


where  ey  is  a  desired  error  tolerance. 

Note  that,  when  j  is  a  large  value,  the  roundoff  errors  due 
to  direct  computations  of  G’  in  eqn.  43  may  occur.  To 
reduce  the  errors,  the  recursive  algorithm  in  eqn.  28a  with 
A  =  (G  — I2n)  (G  +  /2n )_  1  can  be  applied  to  determine  Gj  in 
eqn.  43c  where  Gy  =  sign}0  (A).  Using  the  Gy,  the  approximate 
W  (defined  as  Wj)  yields 


±  Gj  + 1  = 


(Wn)j  (Wn)j 

(IV21 V  <^)/ 
(K^yec"*" 


i,  k  =  1,2 


(44) 


Thus  the  approximate  Riccati  gain  matrix  P  (defined  as 
Pj)  becomes 

^  =  (H'ril/Wi);'  (45a) 

optimal  control  law  and  gain  F  (defined 


(41) 


and  the  approximate 
as  Fj)  become 

u(k)  =  -FjX(k) 

where 

Fj  =  (R  +  BTPjBy'BTPjA 


(45*) 

(45c) 


